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Abstract 

In  this  paper,  we  develop  a  three-dimensional  parallel  solver  using  the  fifth  order  high- 
resolution  weighted  essentially  non-oscillatory  (WENO)  finite  difference  scheme  to  perform 
extensive  simulation  for  three-dimensional  gaseous  detonations.  A  careful  study  is  conducted 
for  the  propagation  modes  of  the  three-dimensional  gaseous  detonation  wave-front  structure 
in  a  long  square  duct  with  different  widths  under  different  initial  perturbations.  The  nu¬ 
merical  results  indicate  that,  with  a  transverse  sinusoidal  perturbation  of  the  initial  ZND 
profile,  when  the  width  of  the  duct  is  less  than  the  cellular  width  (4.5xLi/2),  unstable  det¬ 
onations  can  trigger  a  spinning  motion  in  the  duct.  The  detonation  wave  propagates  in  a 
single-headed  spinning  motion,  with  a  distinctive  “ribbon”  displayed  on  the  four  walls.  In 
this  case,  the  measured  pitch-to-diameter  ratio  is  approximately  3.42,  which  is  slightly  larger 
than  the  theoretically  predicted  value  3.128  for  a  round  duct.  When  the  channel  width  is 
greater  than  the  cellular  width,  detonation  waves  propagate  in  an  out-of-phase  rectangular 
mode.  With  a  transverse  cosine  perturbation  of  the  initial  ZND  profile,  the  front  of  the 
stable  detonation  has  a  rectangular  structure,  and  regular  cellular  patterns  and  in-phase 
“slapping  waves”  can  be  observed  clearly  on  the  four  walls.  The  width-to-length  ratio  of 
the  cellular  patterns  is  approximately  0.5.  For  a  mildly  unstable  detonation,  its  front  has 
an  in-phase  rectangular  structure  at  the  early  stage,  then  the  wave- front  becomes  fiat.  Over 
time,  but  it  still  maintains  an  in-phase  rectangular  structure  after  reigniting.  For  highly 
unstable  detonations,  the  wave-front  has  a  rectangular  structure  at  the  early  stage.  After 
a  low  pressure  stage  for  a  very  long  time,  detonation  occurs  once  again.  At  this  time,  the 
detonation  front  structure  becomes  very  twisted,  and  the  triple-lines  become  asymmetrical. 
Finally,  a  spinning  detonation  mode  is  triggered.  With  a  symmetrical  perturbation  mode 
along  the  diagonals  of  the  detonation  front,  for  the  stable  detonation,  an  diagonal  detonation 
is  formed  and  the  detonation  front  maintains  a  diagonal  structure,  but  no  “slapping  waves” 
appears  on  the  walls.  The  width-to-length  ratio  of  the  cellular  structure  is  equal  to  that  in 
the  rectangular  structure.  For  mildly  unstable  and  highly  unstable  detonations,  the  front 
has  a  diagonal  structure  at  the  early  stage.  After  a  short  period  of  time,  the  diagonal  struc¬ 
ture  of  the  detonation  front  cannot  be  maintained,  and  it  ultimately  evolves  into  a  spinning 
detonation. 
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1  Introduction 


Detonation  is  a  complex  supersonic  flow  phenomenon  where  its  front  consists  of  a  precursor 
shock  wave  that  propagates  into  the  unreacted  medium  at  supersonic  speed  with  a  thin 
reaction  zone  immediately  behind  the  shock.  The  precursor  shock  compresses  the  unreacted 
medium  and  increases  its  temperature.  Burning  occurs  behind  the  front,  which  can  release 
a  large  amount  of  heat  that  in  turn  supports  the  precursor  shock  wave  to  keep  propagating 
forward. 

Numerous  experimental  and  numerical  studies  have  been  performed  to  study  detonation. 
However,  detailed  numerical  studies  on  detonation  mostly  remain  in  two-dimensional  sim¬ 
ulations.  As  is  well  known,  detonation  is  essentially  a  three-dimensional  phenomenon,  and 
some  important  structural  features  cannot  be  obtained  from  two-dimensional  simulations, 
such  as  the  slapping  waves  and  spinning  structures.  Unfortunately,  a  highly  rehned  grid 
resolution  is  required  for  such  computation,  especially  for  unstable  detonation,  which  makes 
three-dimensional  numerical  simulation  of  detonation  tremendously  computing-resource  in¬ 
tensive.  In  this  paper,  we  attempt  to  address  this  difficulty  by  developing  a  high  order 
accurate  weighted  essentially  non-oscillatory  (WENO)  solver  with  parallel  implementation 
so  that  the  desired  wave  structures  can  be  resolved  within  acceptable  computational  time. 
WENO  schemes  have  the  advantage  of  high  order  accuracy  and  robust,  sharp,  and  essentially 
non-oscillatory  shock  resolution  [1,  2,  3,  4]. 

Previous  research  in  the  literature  indicates  that  the  heat  of  the  reaction,  the  activation 
energy,  the  overdrive  factor  and  the  ratio  of  specific  heat  all  have  some  effects  on  the  stability 
of  detonation  and  cellular  patterns  [5,  6].  Bourlioux  et  ah  [7],  Papalexandris  et  ah  [8],  He 
and  Karagozian  [9],  and  Daimon  and  Matsuo  [10]  carried  out  numerical  simulation  of  one¬ 
dimensional  detonations.  Guirguis  et  ah  [11],  Bourlioux  and  Majda  [12],  Papalexandris 
et  al.  [13],  Gamezo  et  al.  [14],  Hwang  et  al.  [15],  and  Shepherd  et  ah  [16]  carried  out 
numerical  simulation  of  two-dimensional  detonations.  These  numerical  simulations  have 
investigated  the  instability  of  the  detonation  in  details  and  have  drawn  some  very  important 
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conclusions.  The  instability  of  the  detonation  wave  leads  to  irregular  cellular  structure 
which  is  caused  by  transversely  developing  of  triple  points  (triple  lines  in  three-dimensional 
detonation  waves)  along  the  front.  For  unstable  detonations,  a  very  small  initial  perturbation 
will  quickly  develop  transversely  with  time.  Therefore,  with  different  reactive  parameters, 
such  as  the  heat  of  reaction,  the  activation  energy  and  the  overdrive  factor,  and  with  a  given 
initial  perturbation  and  geometrical  configuration,  the  detonation  front  structure  may  have 
different  evolution  processes,  and  hence  the  detonation  modes  may  be  different. 

Detailed  experimental  study  on  three-dimensional  gas  detonation  structure  has  been  con¬ 
ducted  as  early  as  in  the  1960s  (White  and  Cary  [17]  and  Strehlow  [18]).  Recent  numerical 
simulations  and  experimental  results  have  shown  that  there  are  three  main  types  of  cellular 
detonation  structures,  namely  rectangular,  diagonal,  and  spinning  modes  [19,  20].  Hanana 
et  al.  [21],  by  experimental  study,  pointed  out  that  in  rectangular  tubes  at  least  two  types 
of  detonation  structure  exist,  namely  the  rectangular  and  diagonal  structures.  Through  the 
soot-foil  tracks  on  the  walls,  they  showed  the  diamond-shaped  cellular  patterns.  They  also 
showed  that  the  rectangular  structure  is  characterized  by  straight  triple  lines  emanating  from 
the  leading  front  that  are  parallel  or  orthogonal  to  the  walls  of  the  flow  domain.  Another 
characteristic  feature  is  the  occurrence  of  slapping  waves  on  the  walls.  These  waves  are 
formed  by  collisions  between  a  triple  line  and  the  wall  or  between  two  triple  lines.  They 
divided  rectangular  detonation  modes  into  in-phase  and  out-of-phase.  When  a  cluster  of 
parallel  triple  lines  collide  with  the  walls  simultaneously,  an  in-phase  rectangular  detonation 
is  formed.  On  the  contrary,  it  will  be  an  out-of-phase  rectangular  detonation.  When  the 
movement  direction  of  the  triple-line  is  along  the  diagonal  line,  the  detonation  front  has  a 
diagonal  structure,  and  the  direction  of  the  transverse  wave  propagation  has  a  45°  angle 
with  the  wall,  while  on  the  wall  the  slapping  waves  would  disappear.  Further,  they  pointed 
out  that,  compared  with  a  rectangular  detonation,  a  diagonal  detonation  has  a  high  front 
pressure,  a  large  average  velocity  and  a  reduced  cell  length,  and  that  the  rectangular  det¬ 
onation  formed  by  the  overlapping  of  two  two-dimensional  detonations.  However,  in  their 
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experimental  research,  they  did  not  observe  out-of-phase  detonation.  They  believed  the  ig¬ 
nition  method  was  the  key  factor  for  the  formation  of  different  detonation  modes.  However, 
they  did  not  give  the  relationship  between  different  detonation  modes,  the  front  features 
of  different  detonation  modes  and  the  mechanism  of  their  generation.  Williams  et  ah  [22] 
adopted  a  one-step  reaction  model  and  performed  numerical  simulation  of  three-dimensional 
detonation  in  a  rectangular  duct,  and  observed  the  front  rectangular  structure.  Tsuboi  et 
ah  [23]  carried  out  three-dimensional  numerical  simulation  of  detonation  propagation  in 
the  rectangular  duct,  and  further  conhrmed  the  detonation  is  divided  into  in-phase  rectan¬ 
gular  and  out-of-phase  modes.  Deledicque  and  Papalexandris  [24]  studied  the  rectangular 
and  diagonal  structures  of  detonations  in  a  three-dimensional  rectangular  tube  by  using  a 
one-step  chemical  reaction  model.  When  the  overdrive  factor  /=  1.1,  the  heat  of  reaction 
Q  =  2.0,  and  the  activation  energy  Ea  =  20,  the  detonation  is  in  a  rectangular  out-of-phase 
mode  under  the  transverse  sinusoidal  perturbation,  and  out-of-phase  “slapping  waves”  are 
formed  on  the  walls  which  are  mutually  perpendicular  in  the  square  tube.  When  they  used 
a  constant  perturbation  along  the  diagonals  of  the  front,  a  diagonal  detonation  was  formed. 
They  also  studied  the  rectangular  mode  of  the  unstable  detonation  for  the  activation  energy 
Ea  =  50,  the  heat  of  reaction  Q  =  50,  and  the  overdrive  factors  /=  1.2  and  /=  1.6  under  a 
sinusoidal  perturbation,  but  they  did  not  study  the  evolution  of  the  two  unstable  detonation 
fronts  in  perturbation  along  the  diagonals  of  the  front. 

Campbell  and  Woodhead  [25]  identihed  the  phenomenon  of  spinning  detonations  in  small- 
diameter  tubes  near  the  detonation  limits.  They  found  that  the  pitch  of  the  spin  is  about 
three  times  the  diameter  of  the  pipe,  see  also  [26].  Lee  et  al.  [27]  observed  a  single-headed 
spinning  detonation  in  a  square  channel  through  the  soot-foil  record,  and  the  single  helical 
trajectory  from  the  four  side  walls  of  the  square  channel  is  unfolded  to  show  the  same  soot 
characteristics  as  in  the  round  duct.  Recently,  studies  have  revealed  more  information  on 
spinning  detonation  in  rectangular  channels  and  round  tubes  [23,  28,  29,  30,  31,  32,  33,  34,  35, 
36].  Zhang  et  al.  [26,  30]  conducted  experimental  studies  on  the  two-phase  flow  detonation 
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in  a  round  tube,  and  showed  that  in  a  stable  propagation  of  detonation,  the  transverse  wave 
has  played  a  leading  role.  Huang  et  ah  [31]  conducted  experiments  on  spinning  detonations 
with  a  detailed  analysis  of  the  shock  structure.  Experimental  results  indicate  that  the  actual 
structure  of  the  spinning  detonation  tries  to  match  closely  to  the  condition  where  the  state 
parameters  (pressure  and  temperature)  reach  their  maximum  values.  Tsuboi  and  co-workers 
[23,  36,  37,  38]  also  carried  out  extensive  numerical  simulations  on  the  detonation  in  round 
tubes  and  rectangular  ducts.  They  pointed  out  that  the  formation  of  an  unburned  gas  pocket 
behind  the  detonation  front  was  not  observed  in  their  results  because  the  rotating  transverse 
detonation  completely  consumed  the  unburned  gas.  Don  et  al.  [39,  40]  investigated  spinning 
detonation  in  narrow  channels  and  detonation  structures  under  different  initial  perturbations 
by  numerical  simulation.  Their  numerical  results  showed  that  the  spinning  detonation  only 
exists  in  narrow  channels.  When  the  channel  width  is  large  enough,  the  spinning  detonation 
goes  away.  These  studies  on  the  spinning  detonation  do  not  consider  modes  of  propagation  of 
detonation  waves  in  tubes  of  different  widths  or  when  different  chemical  reaction  parameters 
are  selected.  Therefore,  with  different  chemical  reaction  parameters,  the  dynamical  behavior 
and  spinning  mechanism  of  the  spinning  detonation  remain  unclear,  and  whether  a  change  in 
duct  width  and  initial  perturbation  can  trigger  spinning  detonation  warrants  further  study. 

The  main  objective  of  this  paper  is  to  simulate  numerically  the  detonation  modes  in 
different  square  tubes,  with  different  chemical  reaction  parameters  and  under  different  initial 
disturbances,  to  provide  a  detailed  front  structure  description  for  each  detonation  mode  and 
further  to  hnd  out  the  conditions  for  triggering  spinning  detonation. 

2  Governing  equations 

The  governing  equations  are  the  three-dimensional  Euler  equations  with  a  source  term  that 
represents  chemical  reactions.  In  conservation  form,  these  equations  may  be  written  in  the 
compact  form 

^  dF{U)  dGjir)  dH{U} 

dt  dx  dy  dz 
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where  the  conserved  variable  vector  U,  the  flux  vectors  F,  G,  and  H  as  well  as  the  source 


term  S  are  given,  respectively,  by 

U  =  (p,  pu,  pv,  pw,  pE,  p  Yf  (2) 

F{U)  =  (pu,  pu^  +  p,  puv,  puw,  pu{E  +  p/ p),  puY)'^  (3) 

G{U)  =  {pv,  pvu,  pi?  +  p,  pvw,  pv{E+  p/ p),  pvY)'^  (4) 

H{U)  =  {pw,  pwu,  pwv,  pv?  +  p,  pw{E+  p/ p),  pwY)'^  (5) 

S{U}  =  {0,0,0,0,0,pujf  (6) 

RT  1 

E= - -+YQ+-{u^  +  v^  +  w^)  (7) 

7  —  i  2 

u  = -KpYe-^^^^^^  (8) 

p  =  pRT  (9) 


here  u,  v,  w  are  the  Cartesian  components  of  the  fluid  velocity  in  the  x,  y,  z  directions, 
respectively,  p  is  the  density,  p  is  the  pressure,  E  is  the  total  energy  per  unit  volume,  T  is 
the  temperature,  and  Y  is  the  reactant  mass  fraction.  Q  is  the  heat  of  reaction,  Ea  is  the 
activation  energy,  7  =  1.2  is  the  specific  heat  ratio,  and  K  is  the  pre-exponential  factor. 

Zeldovich,  von  Neumann  and  Doring  [41]  sought  the  one-dimensional  steady  ZND  an¬ 
alytical  solution  in  the  1940s.  When  we  specify  the  steady  state  solution  in  term  of  the 
dimensionless  primitive  variables  in  unreacted  zone  as  p  =  1,P  =  l,u  =  —D,  the  ZND 
analytical  solution  can  be  given  by  the  formula 


D^{j-1) 

(10) 

(11) 

P  =  ^  +  Di{Y)  +  +  1 

7-1 

(12) 

where 

5(1')  =  ?i(D^  +  1)V(7  +  !)£>=  -  +  pYj)  +  (13) 
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and  the  detonation  velocity  =  {Dcj)‘^f,  with  the  overdrive  factor  /  >  1,  and  Dqj  being 
the  C-J  detonation  velocity,  given  as  follows  [20] 

Dcj  =  -  1)  Q/2  +  -  1)  Q/2  +  7Po/po-  (14) 

3  The  numerical  method 


In  this  section  we  briefly  describe  the  general  framework  of  the  the  (2r-l)-th  order  local  char¬ 
acteristics  based  weighted  essentially  non-oscillatory  (WENO)  conservative  hnite  difference 
scheme  for  solving  the  system  of  hyperbolic  conservation  laws  in  Section  2.  More  details  can 
be  found  in  [2,  3].  In  the  subsequent  computation  we  use  the  hfth  order  (r  =  3)  version. 

Consider  a  uniform  grid  dehned  by  the  grid  points  Xi  =  iAx,  i  =  0,  -  ■  ■  ,71^,  yj  =  jAy, 
j  =  0,  •  •  •  ,ny,  and  Zk  =  kAz,  k  =  0,  •  •  •  where  for  simplicity  Ax  =  Ay  =  is  the 
uniform  grid  spacing.  The  cell  boundaries  are  given  by  XiJ^i/2  =  Xi+Ax/2,  yj+i/2  =  yj+Ay/2, 
Zk+i/2  =  Zk  +  Az/2.  The  original  partial  differential  equation  is  written,  at  the  grid  points. 


in  the  following  form 

dU{xi,yj,Zk,t) 

dt 


dF{U)  dG{U)  dH{U) 


\x=Xi,y=yj,z=Zk 


A  conservative  hnite  difference  scheme  for  approximating  the  equation  (15)  can  be  written 


i— 1/2, j,k  — l/2,fc  f?j,j-|-l/2,A:  Zl^jk— 1/2  -t^i,j,k-i-l/2 


Fiik — 1/2 


+  Si^j^k  (16) 


where  Ui^^k  is  an  approximation  to  the  point  value  U (xj,  yj,  Zk,  t)  of  the  solution.  The  numer¬ 
ical  huxes  Ei±i/2,yfc,  Gij±i/2,k  and  H/j^k±i/2  can  be  computed  using  the  known  grid  values 
of  Fij^k,  Gij^k  and  by  a  WENO  reconstruction.  The  classical  (2r-l)-th  order  WENO 
scheme  uses  a  (2r-l)-point  stencil  to  approximate  the  numerical  hux,  which  is  divided  into 
r  substencils  Sq,  Si, ,  Sr-i  with  each  substencil  containing  r  grid  points. 

Take  as  an  example.  The  WENO  reconstruction  is  made  as  follows.  We  hrst 

perform  a  Lax-Friedrichs  hux  splitting  for  the  hux  T) j  ^  on  the  node,  obtaining 

^tj,k  ~  9  {^{Ui,j,k)  +  OiUij^k)  ,  ^i,j,k  ~  9  {F{Uij^k)  ~  OiUij^k)  (17) 
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where  a  =  maxi<j<„^  \\{Uij^k)\,  for  fixed  j  and  k,  with  \{Uij^k)  being  the  largest  (in  absolute 
value)  eigenvalue  of  the  Jacobian  A{U)  =  It  is  also  possible  to  take  the  maximum  locally 
and/or  to  take  the  maximum  over  individual  eigenvalues,  to  reduce  numerical  dissipation, 
see  [2]  for  more  details.  We  then  construct  also  two  numerical  fluxes  .fo+1/2  jfc  ^i+i/2jk^ 
approximating  and  F~  respectively,  and  their  sum  will  be  the  numerical  flux 
Below  we  will  give  the  details  of  the  construction  of  the  numerical  flux  ^  only,  as  the 

procedure  to  construct  the  numerical  flux  F~-  ^  is  symmetric  to  that  for  FA  ^  with  respect  to 
the  location  i  +  1/2.  The  (2r-l)-th  order  numerical  flux  is  built  through  a  convex 

combination  of  the  r-th  order  numerical  fluxes  F^’W  . .  based  on  the  substencils 

2+1/2, J,/C 

r— 1 


2+1/2, ji,/c  /  ^ 


=  > 


-JO 


2  +  1/2, Jl,/c 


(18) 


z=o 


where  ui  are  the  nonlinear  weights,  satisfying  a;;  >  0  and  =  1-  If  no  discontinuity 

exists  on  the  large  stencil,  there  are  positive  linear  weights  di  which  yield 

r— 1 


2  +  1/2, 


,(0 

+1/2, j,k 


=  F{xi+i/2,  Vj,  Zk)  +  0{Ax  ) 


(19) 


z=o 


where  F{xi+i/2,  Uj,  z^)  is  the  theoretically  infinite  order  accurate  flux.  We  again  refer  to  [2] 
for  more  details. 

By  consistency  =  1.  In  smooth  regions,  we  would  require  oji  =  di  +  0(Aa;’’“^). 

In  WENO  schemes,  the  nonlinear  weight  function  should  be  smooth.  Moreover,  when  the 
substencil  contains  discontinuities,  the  corresponding  weight  ui  should  be  very  small.  The 
weights  in  [2]  are  defined  as 


Ui  = 


ai 


/=  0,...,r-  1; 


d, 


ElJa.’  ■  (■^  +  W 

The  smoothness  indicators  are  given  by 

i-1/2 


(20) 


r^i+i 

P,=Y/  / 

n=l  d  Xi_^/ 


dx^- 


dx 


(21) 


where  FAk(^x,y,  z)  is  the  reconstruction  polynomial  for  the  l-th  substencil,  the  parameter 
e  is  used  to  avoid  the  division  by  zero  in  the  denominator,  and  is  taken  as  10“®  in  the 


simulation.  In  the  case  of  r  =  3,  the  resulting  WENO  scheme  is  hfth-order  accurate,  with 
a  5-point  big  stencil  S^,  and  three  substencils  Sq,  81,82,  each  containing  3-points.  For  the 
substencil  S'o  one  has 


A+,(0)  ^  _rp+ _ p+  _| _ p+ 

i+ll2,j,k  q  i—2,j,k  i—l,j,k  i,j,k 


for  the  substencil  81  one  has 


p+’W  =_1f+  —-F~^  -I--F+  123') 

'^i+l/2,j,k  Qi,j,k'‘^^i+l,j,k 

and  for  the  substencil  82  one  has 

_p+'(2)  =1f+  +-F+  —-F~^  124') 

i+l/2,j,k  2"^  i,j,k  '  Q  Q  *+2,i,fc 

The  linear  weights  are  do  =  3/10,  di  =  3/5  and  ^2  =  1/10.  The  three  smoothness  indicators 


13  1 

A  =  -  -iFUi*  +  FJ‘  +  l(PU,.k  -  -iPUiM  +  3F:,F  (25) 

13  1 

A  =  -  2t;+ ,  +  +  ^{K-i,,k  +  K^+i,j,k?  (26) 

13  1 

-  2E;:i,,,,  +  +  4(3^^+,  -  (27) 

Using  the  same  procedure,  the  numerical  fluxes  in  y  and  2:  directions  k  ^^*4 

,  1  /o  can  also  be  constructed.  We  take 

7/i+l/2J,fc  =  Pj^l/2, j,k  +  ^i+ll2,j,k  (23) 


Gij+ii2,k  -  G'*ti+i/2,fc  +  ^i,j+i/2,k  (29) 

77ij,fc+i/2  =  (30) 

Then  the  numerical  scheme  (16)  is  formed.  If  we  denote  its  right  hand  side  as  L{Uij^k), 

then  the  fully  discrete  scheme,  using  the  third-order  TVD  Runge-Kutta  scheme  [42]  for 

temporal  discretization,  is  given  as 

f  dh  =  ,i) 

{  dh  =  + Aii(dh) 

[  up+l  =  +  AtLiuflP) 
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When  studying  three-dimensional  detonation  problems,  due  to  the  very  large  size  of  the 
computational  domain  and  a  large  number  of  grid  points,  it  is  essential  to  implement  the 
WENO  scheme  in  an  efficient  three-dimensional  parallel  code  to  be  used  on  a  massively 
parallel  platform,  which  is  carried  out  in  this  paper. 

4  Results  and  discussion 

4.1  Grid  convergence  study  on  a  one-dimensional  mildly  unstable 
detonation 

For  this  study  we  take  the  parameters  as  /=  1.6,  Ea  =  50.0,  Q  =  50.0,  and  K  =  230.75 [43]. 
The  size  of  the  computational  domain  is  lOOOxLi/2,  where  T1/2  is  the  length  of  a  half  reaction 
zone.  The  one-dimensional  ZND  analytical  solution  is  used  as  the  initial  condition.  The  left 
boundary  is  an  inflow,  and  the  right  boundary  is  an  outflow.  Three  resolutions  employed  in 
the  simulations  are  5pts/Ti/2  (i-e.  5  points  in  the  length  of  a  half  reaction  zone),  lOpts/ L1/2 
and  20pts j Lij2j  respectively.  The  hnal  time  is  t  =  100.  For  a  one-dimensional  unstable 
detonation,  the  round-off  error  is  enough  to  provide  an  initial  perturbation,  which  further 
triggers  the  physical  instability  of  the  detonation.  Therefore,  it  is  unnecessary  to  artihcially 
add  any  initial  perturbation.  Figure  1  illustrates  the  ZND  front  pressure  time  histories  for  the 
three  grid  resolutions.  It  can  be  seen  from  Figure  1  (a)  that,  when  the  grids  become  hner,  the 
front  pressure  peak-time  curve  tends  to  be  consistent,  and  the  pressure  curves  corresponding 
to  lOpts/Li/2  and  20pts/Li/2  are  largely  overlapping,  indicating  the  numerical  computation 
is  convergent.  Furthermore,  the  advantages  of  the  hfth  order  WENO  scheme  can  be  found, 
i.e.,  the  ZND  pressure  peak  for  lOpts/Li/2  is  98.58;  for  20pts/ L1/2  h  is  98.6,  which  is  in  close 
agreement  with  the  theoretical  value  Pmax  =  98.6  [44].  Comparing  the  curves  in  Figure  1 
(b)  it  can  be  seen  that,  for  the  one-dimensional  mildly  unstable  detonation,  to  reach  the 
theoretical  value  98.6,  it  is  enough  for  the  hfth  order  WENO  scheme  to  use  only  lOpts/Ti/2- 
The  number  of  grids  used  is  obviously  decreased  compared  with  other  algorithms,  e.g.  [45]. 
Moreover,  with  an  increasing  grid  point  number  the  desired  resolution  is  quickly  obtained. 
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a) 


b) 


Figure  1:  Grid  convergence  study  with  the  mildly  unstable  detonation  (overdrive  factor 
/=  1.6).  The  theoretically  predicted  value  is  Pmax  =  98.6. 
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4.2  Study  on  two-dimensional  detonations 

The  propagation  mode  and  cellular  structure  of  two-dimensional  unstable  detonations  with 
different  overdrive  factors  and  different  computational  domain  widths  are  studied  to  provide 
reference  to  the  simulation  of  three-dimensional  detonation  structures.  We  adopt  the  one¬ 
dimensional  ZND  analytical  solution  as  the  initial  condition,  and  add  a  very  small  transverse 
perturbation  to  the  ZND  prohle.  The  left  boundary  is  an  inflow  boundary,  and  the  right 
boundary  is  an  outflow  boundary.  The  upper  and  bottom  boundaries  are  solid  walls. 

4.2.1  Influence  of  different  overdrive  factors  on  two-dimensional  detonation  cell 
widths 

In  order  to  study  the  detonation  cellular  width  with  high  activation  energy  at  different 
overdrive  factors,  the  computational  domain  is  set  as  [240,  4]xLi/2.  The  chemical  reaction 
parameters  are  listed  in  Table  1.  Figure  2  shows  the  maximum  pressure  history  of  the 
detonation  at  different  overdrive  factors.  Due  to  the  very  small  width-length  ratio  of  the 
duct,  the  complete  computational  domain  is  divided  into  three  sections  for  showing  the 
maximum  pressure  history.  It  can  be  seen  from  the  flgure  that,  with  a  decreasing  overdrive 
factor,  the  cell  width  reduces.  When  /  =  1.0,  the  cell  width  is  about  3xLi/2,  the  width- 
length  ratio  of  the  cell  is  0.50.  When  /  =  1.2  and  /  =  1.6,  no  complete  cell  can  be  formed 
because  the  transverse  wavelength  is  larger  than  the  duct  width.  Thus  it  can  be  seen  that, 
with  an  increasing  overdrive  factor,  the  transverse  wavelength  increases. 

Table  1:  Overdrive  factor  and  pre-exponential  factor  [Ea  =  50.0,  Q  =  50.0)  [20,  43] 


A 

B 

C 

Overdrive  factor  / 

1.0 

1.2 

1.6 

Pre-exponential  factor  K 

2566.42 

871.42 

230.75 
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Figure  2:  Maximum  pressure  history  of  detonation  at  different  overdrive  factors,  a)  /  =  1.0, 
b)  /  =  1.2,  c)  /  =  1.6. 

4.2.2  Influence  of  the  two-dimensional  computational  domain  width  on  unstable 
detonation  propagation 

One-dimensional  linear  stability  analysis  suggests  that,  the  detonation  stability  is,  to  a  large 
extent,  dependent  on  the  selected  reaction  parameters.  Lee  and  Stewart  [43]  pointed  out 
that,  at  given  /,  Q  and  7,  with  an  increasing  activation  energy,  the  detonation  becomes 
unstable.  In  this  section,  the  influence  of  the  duct  width  on  the  detonation  propagation 
mode  is  numerically  investigated.  In  the  computation,  we  take  Ea  =  50.0,  Q  =  50.0,  and 
/  =  1.0  as  the  parameters.  The  duct  width  is  taken  as  lxLi/2,  2xLi/2,  4xLi/2,  8xLi/2 
and  16xLi/2  respectively,  and  the  duct  length  is  240xLi/2-  Figure  3  displays  the  maximum 
pressure  history  for  the  computational  domain  of  different  widths.  When  the  width  is  lxLi/2 
and  2xLi/2,  no  complete  detonation  cells  can  be  formed.  When  the  width  is  8xLi/2  and 
16xLi/2,  detonation  cells  can  be  formed,  but  they  are  uneven.  Average  width-length  ratio 
of  the  cells  measured  is  about  0.51,  which  is  slightly  smaller  than  0.55  that  was  measured  in 
[12]. 
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a) 


b) 


d) 


e) 


Figure  3:  Maximum  pressure  history:  a)  lxLi/2,  b)  2xLi/2,  c)  4xLi/2,  d)  8xLi/2,  e) 
16xLi/2. 
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4.3  Study  on  three-dimensional  detonations 

We  use  the  one-dimensional  ZND  analytical  solution  as  the  initial  condition,  and  add  a 
very  small  perturbation  in  the  ZND  prohle  to  accelerate  the  unstable  growth  of  the  flow. 
Unless  otherwise  indicated,  the  left  boundary  of  the  duct  is  an  inflow,  the  right  boundary 
is  an  outflow,  and  the  remaining  sides  are  rigid  walls.  The  detonations  propagate  along 
the  x-direction.  Simulations  of  three-dimensional  detonations  demand  high  grid  resolution, 
especially  for  unstable  detonations.  Therefore,  a  grid  resolution  of  20pts/Li/2  is  employed 
to  cover  the  computational  domain. 

4.3.1  Verification  of  the  numerical  resolution  of  three-dimensional  detonations 

We  add  source  terms  (functions  of  x,  y,  z  and  t)  to  the  PDE  system  so  that  the  following 
functions  are  exact  solutions  to  the  modihed  PDE  system: 

p{x,  y,  z,  t)  =  1  +  I  sin(7r(x  +  y  +  z  -  2t)) 

u(x,  y,  z,  t)  =  v(x,  y,  z,  t)  =  w{x,  y,  z,  t)  =  p{x,  y,  z,t)  =  1 

Y(x,  y,  z,  t)  =  ^  ^  sin(7r(x-  2y+2z+  t))  . 

The  initial  condition  is  set  to  be 

p{x,  y,  z,  t)  =  1  +  I  sin(7r(a:  +  y  +  z)) 
u{x,  y,  z,  0)  =  v{x,  y,  z,  0)  =  w{x,  y,  z,  0)  =  p{x,  y,  z,  0)  =  1 
Y{x,  y,  z,  0)  =  ^  +  j  sin(7r(a:  -2y+  2z))  . 

The  parameter  values  Ea  =  50.0,  /=  1.6,  and  Q  =  50.0  are  given,  and  the  computational 

domain  is  taken  as  [0,  2]x[0,  2]  x  [0,  2].  We  use  uniform  meshes  with  periodic  boundary 

conditions,  and  give  the  errors  and  orders  of  accuracy  for  the  density  at  time  t  =  2.0  in  Table 
2.  It  can  be  observed  clearly  that  the  designed  hfth  order  accuracy  is  achieved. 
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Table  2:  and  errors  and  orders  of  accuracy  for  the  density 


h 

CFL 

^00  Error 

Order 

Error 

Order 

Error 

Order 

2/20 

0.2 

1.26x10-^ 

8.13x10-^ 

8.96x10-® 

2/40 

0.2 

4.48x10-^ 

4.84 

2.46x10-® 

4.93 

2.74x10-® 

5.05 

2/80 

0.2 

1.32x10-® 

5.23 

7.50x10-^ 

5.05 

8.12x10-^ 

5.09 

2/160 

0.2 

4.35x10-® 

4.95 

2.15x10-® 

5.14 

2.44x10-® 

5.08 

2/320 

0.2 

1.16x10-9 

5.25 

6.92x10-1° 

4.98 

9.01x10-1° 

4.78 

4.3.2  Influence  of  the  duct  width  on  three-dimensional  detonations 

From  the  numerical  results  for  two-dimensional  detonations,  we  have  seen  that,  when  /=  1.0, 
Q  =  50.0,  and  Ea  =  50.0,  the  transverse  wavelength  of  the  unstable  detonation  is  very  small, 
and  the  cellular  width  is  relatively  small.  It  is  appropriate  to  choose  this  set  of  chemical 
reaction  parameters  to  study  the  influence  of  three-dimensional  duct  width  on  detonation 
propagation  modes.  We  study  three-dimensional  square  ducts,  with  length  160xLi/2,  and 
with  width  lxLi/2  (Case  1),  2xLi/2  (Case  2),  4xLi/2  (Case  3),  8xLi/2  (Case  4)  and  16xLi/2 
(Case  5),  respectively.  From  our  experience  in  the  grid  convergence  study  on  one-dimensional 
unstable  detonation,  we  can  see  that  lOpts/Li/2  can  meet  the  required  computational  accu¬ 
racy.  However,  for  three-dimensional  unstable  detonations,  more  number  of  grid  points  is 
required  for  the  half  reaction  zone  [46,  6,  47].  Hence,  20pts/ L1/2  is  adopted  in  our  computa¬ 
tion.  We  use  the  one-dimensional  ZND  analytical  solution  as  the  initial  condition,  and  add 
a  very  small  transverse  sinusoidal  perturbation  to  the  ZND  prohle.  Parallel  computation 
is  performed  by  using  high-performance  computers.  The  computational  domain,  number  of 
grid  points  and  number  of  MPI  processors  are  shown  in  Table  3. 


Table  3:  Computational  domain,  number  of  grid  points  and  number  of  MPI  processors 


Case 

1 

2 

3 

4 

5 

Domain 

[160,1,1]xLi/2 

[160,2,2]  xLi/2 

[160,4,4]  X  Li/2 

[160,8,8]  X  Li/2 

[160,16,16]  xLi/ 

Grids 

3200x20x20 

3200x40x40 

3200x80x80 

3200x160x160 

3200x320x320 

Processors  64 

180 

270 

360 

450 

Figure  4  shows,  in  the  duct  with  width  2xLi/2,  the  typical  front  feature  and  density 
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gradient  on  the  walls  of  the  detonation  in  a  cycle.  The  front  structure  is  similar  to  that  in 
[48,  49],  suggesting  spinning  detonations  have  formed.  Due  to  the  differences  in  the  type 
of  perturbation,  spinning  detonation  rotates  counter-clockwise,  opposite  to  that  in  [48].  It 
can  be  seen  that,  the  triple  lines  and  transverse  waves  collide  with  the  walls,  and  strong 
explosions  take  place  near  the  walls.  Transverse  waves  play  an  important  role  in  changing 
the  spinning  direction. 

^  ^ 

a)  b)  c)  d)  e) 

f)  g)  h)  i)  j) 

Figure  4:  Typical  spinning  detonation  density  contour  of  the  detonation  front  for  a  full  cycle. 
Frames(a)-(j)  are  equally  spaced  with  the  same  time. 

Figure  5  (a,  b,  c)  displays  the  maximum  pressure  history  and  the  main  spinning  tracks  on 
the  walls  when  the  detonations  for  Case  1,  Case  2  and  Case  3  appear  spinning  mode.  It  can 
be  clearly  seen  from  these  pictures  that  the  spinning  tracks  of  the  single-headed  detonations 
on  the  walls  propagate  counter-clockwise  along  the  walls.  For  Case  1,  the  main  spinning 
dominates.  Because  transverse  waves  are  restricted  by  the  walls,  they  are  not  completely 
developed.  Colliding  with  the  duct  corner,  the  triple  lines  are  not  strongly  reflected.  Hence, 
the  tracks  of  the  reflected  triple  lines  are  not  apparent.  After  the  stable  single- headed 
spinning  detonation  is  formed  for  a  time  period,  its  pressure  decreases,  and  ultimately  it 
cannot  stably  propagate.  With  width  increasing,the  spinning  detonation  formed  can  stably 
propagate  and  the  tracks  of  the  reflected  triple  points  become  more  obvious,  as  shown  in 
Figure  5(c).  Figure  6  shows  the  maximum  pressure  history  on  the  unfolded  walls  for  Case 
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1,  Case  2  and  Case  3.  It  can  be  seen  from  this  figure  that,  for  Case  2  and  Case  3,  a  stable 
single-headed  spinning  detonation  is  formed.  The  main  spiral  and  the  tracks  of  reflective 
triple  point  characterized  by  single-headed  spinning  detonations  [36,  27]  are  shown  on  the 
rigid  walls.  Comparing  the  computational  results  for  Case  2  and  Case  3,  we  can  see  that, 
when  the  width  is  doubled,  the  cycle  of  the  spinning  detonation  also  doubles.  [50,  51]  gave 
the  theoretical  value  of  the  pitch-to-diameter  ratio  of  the  spinning  detonation  in  a  round 
duct  as 

^  _  7r(7  -F  1) 

d  ~  1.8417 

where  pi  is  the  pitch  of  spin,  and  d  is  the  duct  diameter.  The  above  formula  indicates  that 
the  pitch-to-diameter  ratio  is  related  only  to  the  specihc  heat  of  gas.  For  the  gas  having 
7=1.2,  the  pitch-to-diameter  ratio  is  3.128.  Although  the  theoretical  value  was  obtained  in 
the  round  duct,  it  should  be  applicable  also  to  a  rectangular  duct.  Through  the  maximum 
pressure  history  on  the  walls,  we  can  measure  the  pitch-to-diameter  ratio  of  Case  3  as  3.42, 
which  is  slightly  larger  than  the  theoretical  value,  and  also  larger  than  2.7,  the  pitch-to- 
diameter  ratio  measured  in  [26].  The  measured  spinning  angle,  which  can  be  computed  by 
the  length  covered  by  the  detonation  front  along  the  x-direction  in  a  cycle  and  the  duct 
perimeter,  is  51.5°  for  Case  2,  and  is  49.5°  for  Case  3.  They  are  in  good  agreement  with 
the  results  in  [49,  52].  Therefore,  width  change  exerts  little  influence  on  the  spinning  angle. 
When  the  width  increases  to  the  one  for  Case  4,  irregular  cells  are  shown  on  the  walls,  as 
shown  in  Figure  7.  When  the  width  is  further  increased  to  Case  5,  transverse  perturbations 
are  little  affected  by  the  width.  Over  time,  obvious  out-of-phase  slapping  waves  are  displayed 
on  the  walls.  Detonation  gradually  evolves  into  an  out-of-phase  rectangular  mode, as  shown 
in  Figure  8.  The  feature  coincides  with  that  in  [24].  The  measured  cellular  width  is  about 
4.5xLi/2,  larger  than  that  in  the  two-dimensional  computational  results  shown  in  Figure  3 
(c).  It  seems  that  4.5xLi/2  is  the  critical  width  for  the  formation  of  spinning  detonations. 
When  the  duct  width  is  smaller  than  this  value,  spinning  detonations  can  be  triggered.  For 
wider  ducts,  the  wall  boundary  conditions  exert  little  influence  on  the  transverse  waves. 
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so  the  measured  cell  width  can  be  believed  as  the  true  width  of  the  unstable  detonation. 


Comparing  the  computational  results  for  Cases  1-5,  we  can  conclude  that,  when  the  duct 
width  is  smaller  than  the  cell  width,  the  unstable  detonation  propagates  in  a  spinning  mode; 
when  the  duct  width  is  larger  than  the  cell  width,  spinning  detonation  disappears. 


c) 

Figure  5:  Maximum  pressure  history  and  main  spinning  tracks  on  the  walls  (Detonation 
wave  propagates  along  the  x-direction):  a)  Case  1,  b)  Case  2,  c)  Case  3. 


In  order  to  study  the  influence  of  the  overdrive  factor  on  the  spinning  detonation,  we 
take  /=  1.2,  Q  =  50.0  and  Ea  =  50.0.  The  duct  width  is  the  same  as  the  one  for  Case  2. 
Figure  9  shows  the  maximum  pressure  history  of  the  highly  unstable  detonation  on  the  duct 
walls.  Comparing  Figure  6  (b)  with  Figure  9,  we  can  see  that  the  two  unstable  detonations 
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a) 


b) 


Figure  6:  Maximum  pressure  history  of  the  spinning  detonation  on  the  walls:  a)  Case  1,  b) 
Case  2,  c)  Case  3  (detonation  wave  propagates  from  left  to  right). 


Figure  7:  Maximum  pressure  history  on  the  walls:  Case  4. 
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Figure  8:  Maximum  pressure  history  on  the  walls:  Case  5. 


Figure  9:  Maximum  pressure  history  of  the  spinning  detonation  on  the  walls: 
Q  =  50.0,  Ea  =  50.0  (detonation  wave  propagates  from  left  to  right). 
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can  trigger  the  spinning  mode  in  the  duct  with  width  2xLi/2,  and  the  measured  spinning 
angle  and  pitch-to-diameter  ratio  in  the  two  cases  are  equal.  This  is  because  the  transverse 
wavelength  at  /  =  1.2  is  larger  than  that  at  /  =  1.0,  which  has  been  obtained  from  the 
two-dimensional  computational  results, as  shown  in  Figure  2. 

4.3.3  Detonation  structure  at  different  chemical  reaction  parameters  under 
different  types  of  perturbation 

With  different  chemical  reaction  parameters,  the  three-dimensional  detonation  front  features 
and  propagation  mode  under  transverse  cosine  and  symmetrical  perturbation  along  the  di¬ 
agonal  direction  is  investigated.  Give  three  groups  of  chemical  reaction  parameters:  (1) 
Ea  =  20.0,  Q  =  2.0,  /  =  1.1  and  K  =  1134363.64,  corresponding  to  the  one-dimensional 
problem  which  is  a  stable  detonation  [43];  (2)  Ea  =  50.0,  Q  =  50.0,  /=  1.6  and  K  =  230.75, 
corresponding  to  the  one-dimensional  problem  which  is  a  mildly  unstable  detonation  [53]; 
(3)  Ea  =  50.0,  Q  =  50.0,  /  =  1.2  and  K  =  871.42,  which  corresponds  to  a  case  with  high 
activation  energy  and  high  heat  of  reaction,  and,  according  to  a  linear  stability  analysis 
in  Erpenbenk  [53]  and  Lee  and  Stewart  [43],  the  corresponding  one-dimensional  problem  is 
unstable.  Erpenbenk  [54]  also  proposed  that  the  corresponding  two-dimensional  problem  is 
a  highly  unstable  detonation  for  this  case  (3). 

The  initial  condition  is  the  one-dimensional  ZND  analytical  solution,  and  transverse 
cosine  perturbation  or  symmetrical  perturbation  along  the  diagonal  direction  is  added  to 
the  ZND  prohle.  We  take  the  computational  domain  in  length,  width  and  height  as  [144, 
9.6,  9.6]xLi/2.  a  grid  with  2880x192x192  points  are  used. 

(1)  The  detonation  structure  under  transverse  cosine  perturbation 

•  Stable  detonation  {Ea  =  20.0,  Q  =  2.0,  /=  1.1) 

We  take  5pts/ Lij2,  lOpts/Li/2  and  20pts/ L1/2  to  verify  the  grid  resolution  and  the  conver¬ 
gence  of  the  numerical  method.  Figure  10  gives  the  pressure  gradient  distribution  at  f  =  60.0 
on  the  wall  (?/  =  0).  It  can  be  seen  that,  the  hner  the  grid,  the  clearer  the  front  structure. 
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At  bpts/Lii2,  the  global  features  of  the  flow  can  be  seen  only  roughly.  At  lOpfs/Li/2,  the 
detonation  front  structure  and  the  flow  features  have  been  reasonably  well  captured,  how¬ 
ever  the  image  resolution  is  not  high.  At  20pts/Z/i/2,  the  detonation  front  structure  and  flow 
features  are  captured  very  well. 

Figure  11  records  the  maximum  pressure  history  on  the  walls.  It  can  be  seen  from  this 
hgure  that  the  triple-line  trajectories  form  regular  cellular  structures.  Detonation  cells  are 
regular.  The  width-length  ratio  of  the  cell  is  measured  to  be  0.51,  which  is  close  to  the  two- 
dimensional  result  measured  in  [12,  13].  The  slapping  waves  orthogonal  to  each  other  are 
displayed  on  the  walls.  The  slapping  waves  are  parallel  to  the  y  and  ;2- direct  ions,  respectively, 
which  is  in  agreement  with  the  result  in  [24].  Figure  12  displays  the  front  structure  and  the 
pressure  gradient  on  the  walls.  It  can  be  seen  in  this  hgure  that  the  detonation  wave  has 
a  rectangular  structure,  and  the  front  contains  four  pairs  of  triple  lines,  two  of  which  are 
parallel  to  the  y-direction  and  the  other  two  parallel  to  the  ;2-direction.  When  each  pair 
of  the  triple  lines  moves  in  opposite  directions  parallel  to  the  front,  a  pair  of  triple  lines 
will  collide  with  the  walls  or  with  each  other.  When  the  triple  line  collides  with  the  walls, 
slapping  waves  are  shown  on  the  walls.  But  when  a  pair  of  triple  lines  collide  with  each  other, 
the  triple  line  parallel  to  the  propagation  direction  will  appear.  Because  the  perturbation 
mode  is  in  phase,  the  triple  lines  parallel  to  the  ^/-direction  and  to  the  ;2-direction  arrive  at 
the  walls  at  the  same  time,  and  they  collide  with  the  walls.  In-phase  slapping  waves  appear 
on  the  walls.  Hence,  under  transverse  cosine  perturbation,  detonation  waves  propagate  in 
in-phase  rectangular  mode.  Figure  13  displays  the  density  gradient  and  front  structure  on 
the  walls.  From  this  figure,  the  triple  line  movement  direction  and  their  collision  with  the 
walls  can  be  observed  more  clearly. 

•  Mildly  unstable  detonation  (Ea  =  50.0,  Q  =  50.0,  /=  1.6) 

We  take  5pts/Lii2,  lOpts/Li/2  and  20pts/Lij2  to  verify  the  grid  resolution  and  conver¬ 
gence  of  the  mildly  unstable  detonation  under  transverse  cosine  perturbation.  Figure  14 
gives  the  pressure  gradient  distribution  on  the  wall  {y  =  0)  at  t  =  40.0.  It  can  be  seen  that 
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c) 

Figure  10:  Grid  resolution:  a)  5pts/Fi/2,  b)  lOpts/Fi/2,  c)  20pts/Z/i/2- 
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Figure  11:  Maximum  pressure  history  on  the  walls. 


a)t=21.5  b)t=23.5  c)t=25.5  d)t=27.5 

Figure  12:  Front  structure  and  pressure  gradient  on  the  walls. 


a)t=20.5  b)t=22.5  c)t=23.5  d)t=25.5 

Figure  13:  Rectangular  front  structure  and  density  gradient  on  the  walls. 
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at  20pts/Li/2,  the  detonation  front  structure  and  flow  features  can  be  captured  very  well. 

Figure  15  displays  the  evolution  of  a  mildly  unstable  detonation  front  and  the  density 
gradient  on  the  walls.  At  t=  1.46,  due  to  flow  instability,  local  overdrive  is  obvious,  leading 
to  a  highly  distorted  detonation  front.  At  t  =  6.87  —  10.62,  detonations  enter  into  a  low 
velocity  stage,  the  particles  penetrating  into  the  front  need  a  very  long  induction  period. 
The  front  becomes  flat,  the  front  thickness  increases,  and  combustion  is  incomplete.  Over 
time,  detonations  re-appear.  Due  to  the  nature  of  such  detonation  instability,  the  triple  lines 
are  relatively  bent,  and  the  front  shows  “convex”  and  “concave”  features.  At  f  =  13.29,  the 
front  is  relatively  distorted,  but  still  maintains  a  rectangular  structure.  The  triple  lines  on 
the  front,  parallel  to  the  x-  or  the  ^/-directions,  move  toward  the  walls,  and  collide  with  the 
walls  and  then  are  reflected,  moving  at  opposite  directions.  At  f  =  23.49,  they  collide  once 
again  on  the  central  line.  For  three-dimensional  detonations,  the  concave  part  on  the  front 
is  equivalent  to  the  Mach  stem  in  the  two-dimensional  case.  The  flow  between  the  concave 
part  and  other  concave  parts  is  incident  waves.  At  t  =  32.11  —  33.48,  the  triple  lines  collide 
with  the  walls.  Local  explosions  occur  at  the  wall  corner,  and  the  front  is  very  distorted, 
leading  to  inconsistent  reaction  rate  behind  the  front.  Hence,  a  large  quantity  of  unreacted 
pockets  appear  behind  the  front  (Figure  16  (a,b)).  It  can  be  seen  from  Figure  16  (b)  that 
in-phase  slapping  waves  appear  on  the  walls,  and  become  thick  and  bent.  This  is  because  of 
the  flow  instability  and  high  overdrive  leading  to  an  increase  of  unreacted  substances.  The 
measured  width-length  ratio  of  the  cells  is  0.46,  which  is  consistent  with  the  result  in  [24]. 

•  Highly  unstable  detonation  {Ea  =  50.0,  Q  =  50.0,  /=  1.2) 

For  three-dimensional  highly  unstable  detonations,  the  precursor  shock  peak  can  reach 
1.5  times  the  stable  ZND  peak,  and  the  length  of  the  reaction  zone  signihcantly  decreases  [51]. 
Therefore,  higher  grid  resolution  is  needed.  We  take  5pts/Li/2,  lOpts/Li/2  and  2t)pts/ Lij2  to 
verify  grid  resolution  and  convergence  of  the  numerical  method.  Figure  17  gives  the  pressure 
gradient  distribution  of  such  detonations  at  t  =  10.0  on  the  walls  {y  =  0).  It  can  be  seen 
that,  as  the  grids  are  rehned,  the  numerical  results  show  good  convergence,  and  the  front 
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c) 

Figure  14:  Grid  resolution:  a)  5pts/Fi/2,  b)  lOpts/Fi/2,  c)  20pts/Z/i/2- 
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e)t  =  13.29  f)t  =  18.84  g)t  =  23.49  h)t  =  23.86 


i)t  =  32.11  j)t  =  32.56  k)t  =  33.02  l)t  =  33.48 


m)t  =  34.06  n)t  =  34.42  o)t  =  34.84  p)t  =  35.22 


Figure  15:  Front  structure  and  density  gradient  on  the  walls  at  different  times. 


Unbum<:  ’ 
Pockets 


a) 


b) 


Figure  16:  Unreacted  pockets  behind  the  detonation  front. 


structure  becomes  clearer.  At  5pts/Li/2,  the  grids  are  too  coarse  to  resolve  global  features 
of  the  flow  for  this  case.  At  lOpts/ L1/2,  the  global  flow  features  of  the  detonations  can 
be  roughly  captured,  but  the  meso  features  of  the  front  cannot  be  clearly  captured.  At 
20pts/ Li/2,  the  front  structure  and  flow  features  of  detonations  can  be  captured  very  well. 
It  seems  that,  for  highly  unstable  detonations,  20pts/Li/2  can  describe  the  front  structure 
well  by  our  numerical  method. 

Figure  18  gives  the  front  structure  of  the  highly  unstable  detonation  and  the  density 
gradient  at  different  times.  It  can  be  seen  that  the  detonation  front  has  a  rectangular  struc¬ 
ture  at  the  early  stage.  At  t  =  10.71,  the  front  becomes  flat  and  thick,  and  the  detonation 
velocity  decreases.  After  a  low  velocity  stage  for  a  long  period  of  time,  at  f  =  24.93  —  28.43, 
detonations  reappear,  at  which  time  the  front  still  has  a  rectangular  structure.  Highly  un¬ 
stable  flow  leads  to  a  big  local  overdrive,  distorted  front  and  asymmetrical  triple  lines.  At 
t  =  39.30  —  42.61,  many  triple  lines  appear  on  the  front,  and  the  front  structure  becomes 
complicated.  With  continuous  evolution  of  the  detonation  front  structure,  the  front  shows 
obvious  features  of  the  spinning  mode,  as  shown  in  Figure  18  at  f  =  52.82  —  54.97.  so,  this 
detonation  can  eventually  trigger  spinning  detonations.  By  the  maximum  pressure  history 
on  the  walls  it  can  be  seen  that,  before  the  spinning  mode  is  formed,  detonation  is  in  the 
rectangular  mode  at  the  early  stage,  and  obvious  slapping  waves  appear  on  the  walls.  Then 
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c) 

Figure  17:  Grid  resolution:  a)  5pts/Fi/2,  b)  lOpts/Fi/2,  c)  20pts/Z/i/2- 
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the  detonation  enters  a  low  velocity  and  low  pressure  stage  for  a  long  time.  At  last,  on  the 
location  at  about  x=  124.8,  detonation  reappears,  and  the  trajectories  of  the  triple  lines  on 
the  walls  become  irregular, as  shown  in  Figure  19  (c). 

(2)  The  detonation  structure  under  symmetrical  perturbation  along  the  diagonal  direction 

•  Stable  detonation  {Ea  =  20.0,  Q  =  2.0,  /=  1.1) 

Figure  20  records  the  maximum  pressure  history  of  detonations  on  the  walls.  It  can  be 
seen  that  the  slapping  waves  on  the  walls  disappear.  This  phenomenon  is  in  agreement  with 
the  result  in  [24].  It  is  worth  noting  that,  along  the  propagation  direction  of  detonations, 
the  cell  size  on  the  walls  increases.  Such  phenomenon  was  not  observed  in  [24].  The  ratio 
between  the  cell  width  on  the  walls  (x  =  28.8  —  48.0)  in  Figure  20  (b)  and  the  one  measured 
in  Figure  12  is  about  0.7,  which  is  in  agreement  with  the  result  in  [24].  Figure  21  (a,b,c,d) 
displays  the  typical  diagonal  structure  of  the  detonation  front,  which  is  in  good  agreement 
with  the  structure  in  [24].  It  can  be  seen  that,  moving  along  the  diagonal  direction,  triple 
lines  do  not  collide  perpendicularly  with  the  walls.  So,  no  slapping  waves  appear  on  the  walls. 
Over  time,  detonation  waves  still  propagate  in  the  diagonal  mode,  but  the  triple  lines  on 
the  front  collide,  coalesce  and  form  the  front  structure  as  shown  in  Figure  21  (e,f).  Because 
the  triple  lines  coalesce  and  the  number  of  triple  lines  moving  along  the  diagonal  direction 
decreases,  the  cycle  for  the  same  pair  of  triple  lines  to  collide  will  increase  in  time.  The 
cell  size  eventually  increases,  as  shown  in  Figure  20  (a).  Therefore,  for  stable  detonations, 
the  detonations  propagate  in  an  in-phase  diagonal  mode,  and  the  front  maintains  a  diagonal 
structure. 

•  Mildly  unstable  detonation  (Ea  =  50.0,  Q  =  50.0,  /=  1.6) 

Figure  22  shows  the  front  structure  at  different  times.  It  can  be  seen  that  the  detonation 
wave  propagates  in  the  diagonal  mode  at  the  early  stage.  At  t  =  32.24,  the  front  becomes  an 
irregular  structure.  Over  time,  the  front  displays  a  structure  featuring  spinning  detonations. 

•  Highly  unstable  detonation  {Ea  =  50.0,  Q  =  50.0,  /=  1.2) 
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a)t  =  2.47  h)t  =  6.68  c)t  =  10.71  d)t  =  14.96 


e)t  =  24.93  f)t  =  25.79  g)t  =  27.65  h)t  =  28.43 


m)t  =  52.82  n)t  =  53.46  o)t  =  54.02  p)t  =  54.97 


Figure  18:  Front  structure  and  density  gradient  on  the  walls  at  different  times. 
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a) 


b) 


c) 


Figure  19:  Maximum  pressure  histories  on  the  wall  |/  =  0:  a)  Ea  =  20.0  ,(5  =  2.0,  /=  1.1, 
b)  Ea  =  50.0,  Q  =  50.0,  /=  1.6,  c)  Ea  =  50.0,  Q  =  50.0,  /=  1.2. 


a) 


b) 


Figure  20:  Maximum  pressure  history  on  the  walls:  b)  is  a  magnihed  picture  of  the  cellular 
structure  indicated  by  the  red  line  in  a). 
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a)t=18.5 


b)t=20.5 


c)t=22.5 


d)t=24.5  e)t=30.5  f)t=34.5 

Figure  21:  Front  structure  and  pressure  gradient  on  the  walls. 


a)t  =  3.82  h)t  =  32.24  c)t  =  53.68  d)t  =  72.71 


Figure  22:  Front  structure  at  different  times. 
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Figure  23  displays  the  evolution  process  of  the  front  structure  at  different  times.  It  can 
be  seen  from  the  hgure  that,  at  f  =  5.08,  the  triple  lines  move  along  the  diagonal  direction, 
and  they  collide  at  the  main  diagonal  (TLl)  and  the  second  main  diagonal  (TL2).  The 
detonation  front  is  very  thin  and  shows  a  diagonal  structure.  The  triple  line  shape  is  shown 
in  Figure  23  (b).  At  t  =  7.80,  the  triple  lines  collide  on  the  central  line  between  TLl  and  TL2. 
The  front  structure  is  shown  in  Figure  23  (e).  At  t=  8.06,  the  front  becomes  flat  and  thick, 
the  reaction  is  incomplete,  and  the  detonation  velocity  decreases.  At  t=  66.83  —  71.34,  the 
front  structure  is  simplihed  as  the  features  shown  in  Figure  24  (c-j).  At  t=  66.83,  the  front 
is  composed  of  three  Mach  stems  and  one  incident  wave.  As  the  incident  waves  advance, 
triple  lines  move  at  the  direction  as  shown  in  Figure  24  (c).  At  t=  67.48,  the  triple  lines  are 
located  on  the  walls.  So,  the  front  is  also  composed  of  two  Mach  stems,  which  is  in  agreement 
with  the  result  in  [36].  The  transverse  waves  collide  with  the  walls  and  are  reflected,  and  the 
triple  lines  move  downward.  At  t=  68.13,  the  front  structure  is  shown  in  Figure  23  (k),  and 
the  triple  lines  move  at  the  direction  as  shown  in  Figure  24  (e).  Hence,  for  highly  unstable 
detonation,  under  symmetrical  initial  perturbation  along  the  diagonal  direction,  detonations 
propagate  in  the  in-phase  diagonal  mode  at  the  early  stage,  and  the  front  has  a  diagonal 
structure.  Over  time,  the  diagonal  structure  of  the  front  breaks  down,  the  front  shows  a 
structure  featuring  spinning  detonations,  and  the  detonations  show  a  transition  from  the 
diagonal  mode  to  the  spinning  mode.  In  Figure  25,  the  measured  spinning  angle  is  50.2°, 
and  the  pitch-to-diameter  ratio  is  3.15.  These  are  very  close  to  the  theoretical  values  [51]and 
numerical  results  [39]. It  is  thus  seen  that,  the  critical  width  for  detonations  at  different 
chemical  reaction  parameters  triggering  spinning  mode  is  different.  For  the  highly  unstable 
detonation,  when  the  duct  width  is  9.6xTi/2,  detonation  can  still  evolve  into  the  spinning 
mode.  For  a  unstable  detonation  {Ea  =  50.0,  Q  =  50.0  and  /=  1.0),  a  spinning  structure 
cannot  be  triggered  in  a  duct  with  width  larger  than  4.5xTi/2  (as  shown  in  Figures  5,  7,  8). 
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a)t  =  4.80  b)t  =  5.08  c)t  =  5.63  d)t  =  6.68 


e)t  =  7.80  f)t  =  8.06  g)t  =  34.89  h)t  =  43.94 


i)t  =  66.83  j)t  =  67.48  k)t  =  68.13  l)t  =  68.76 


m)t  =  69.41  n)t  =  70.05  o)t  =  70.70  p)t  =  71.34 
Figure  23:  Front  structure  at  different  times. 
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Figure  24:  Schematic  of  the  spinning  detonation  front  structure:  TLl  is  the  main  diagonal 
triple  line,  TL2  is  the  second  main  diagonal  triple  line.  I  is  incident  wave,  M  is  “Mach  leg” . 
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b) 


Figure  25:  Maximum  pressure  history  on  the  walls:  a)  y  =  0.0,  h)  z  =  0.0. 

5  Concluding  remarks 

Using  the  high-resolution  hfth  order  WENO  scheme  with  third  order  TVD-Runge-Kutta  tem¬ 
poral  discretization,  we  perform  extensive  numerical  simulation  for  the  propagation  modes 
of  three-dimensional  gaseous  detonation  front  structure  in  long,  straight  and  rectangular 
ducts  with  different  width  and  under  different  initial  perturbations.  The  following  conclu¬ 
sions  are  reached:  (1)  Under  a  transverse  sinusoidal  perturbation,  when  the  duct  width  is 
smaller  than  the  cellular  width  (4.5xLi/2),  the  unstable  detonation  can  trigger  the  spin¬ 
ning  mode  in  the  duct.  The  detonation  wave  propagates  in  a  single-headed  spinning  mode, 
and  the  ribbon  features  of  single-headed  spinning  detonations  are  shown  on  the  walls.  The 
measured  pitch-to-diameter  ratio  is  about  3.42,  which  is  a  little  larger  than  the  theoretical 
value  3.128.  When  the  duct  width  is  larger  than  the  cellular  width,  the  detonation  tran¬ 
sitions  to  an  out-of-phase  rectangular  mode.  (2)  Under  a  transverse  cosine  perturbation, 
the  stable  detonation  front  has  a  rectangular  structure,  obvious  slapping  waves  appear  on 
the  walls,  and  the  measured  cellular  width-length  ratio  is  about  0.50.  The  mildly  unstable 
detonation  front  has  a  rectangular  structure  at  the  early  stage.  Over  time,  the  front  becomes 
flat,  and  at  last  it  becomes  a  rectangular  mode.  The  highly  unstable  detonation  front  has  a 
rectangular  structure  at  the  early  stage,  but  after  a  low  velocity  stage  for  a  long  period  of 
time,  it  reignites  and  forms  detonation,  the  detonation  front  becomes  very  distorted  and  the 
detonation  becomes  unstable.  Eventually,  the  spinning  detonation  is  triggered.  (3)  Under 
the  symmetrical  perturbation  along  the  diagonal  of  the  front,  the  stable  detonation  is  in  an 
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in-phase  diagonal  mode,  the  detonation  front  maintains  a  diagonal  structure,  the  slapping 
waves  on  the  walls  disappear,  and  the  measured  cell  width-length  ratio  is  equal  to  that  in  the 
case  of  rectangular  structures.  Mildly  unstable  detonation  and  highly  unstable  detonation 
fronts  have  a  diagonal  structure  at  the  early  stage,  but  such  structure  is  very  unstable.  After 
a  short  period  of  time,  the  diagonal  structure  of  the  front  breaks  down  and  eventually  it 
evolves  into  a  spinning  detonation.  Therefore,  spinning  detonation  is  the  ultimate  mode  of 
detonation.  Triggering  spinning  detonation  is  of  significance  for  stable  propagation  of  highly 
unstable  detonations. 
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